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ABSTRACT 

In this paper, we modify Laskar's simplified model of planetary evolution and accre- 
tion to account for the full conservation of the total angular momentum of the system, 
and extend it to incorporate an accretion probability that depends on the mass and 
relative velocity of the colliding particles. We present statistical results for the mass 
and eccentricity of the planets formed, in terms of their semi-major axes, for a large 
number of realizations of different versions of the model. In particular, we find that by 
combining the mass-dependent accretion probability and the velocity-selection mech- 
anism, the planets formed display a systematic occurrence at specific locations. By 
introducing properly scaled variables, our results are universal with respect to the 
total angular momentum of the system, the mass of the planetesimal disc, and the 
mass of the central star. 

Key words: celestial mechanics - planets and satellites: formation - methods: nu- 
merical - methods: statistical 



1 INTRODUCTION 

Essentially, all models are wrong, but some are useful. 

Box fc Draper (1987) 

Recent advances in the understanding of the final 
stages of the formation of the Solar System within the 
Nice model, have achieved, among others, reproducing 
the a rchitecture of the outer part to an unprecede nted 
level l|Tsiganis et al. I l2005l: iMorbidelli et al.ll2007l . 120091 ). In 
the 'first' Nice model |Tsiganis et al.ll2005l ), the initial con- 
ditions are defined such that Jupiter is close to its ac- 
tual position, Saturn is placed just below the semi-major 
axis corresponding to the 1:2 mean-motion resonance with 
Jupiter, and the icy giants are placed in the intervals 11- 
13 AU (Uranus) and 13.5-17 AU (Neptune); the planetary 
masses are the current ones. In addition, the Nice model 
assumes a 30-50 Mq planetesimal disc placed beyond the 
initial orbit of Neptune which ends between 30-35 AU. The 
planetesimal disc transfers angular momentum to the plan- 
ets m aking them migrate ijFernandez fc id 11984 : iMalhotral 
1 19951 ). allowing Jupiter and Saturn to cross their mutual 
1:2 mean-motion resonance. This event increases their ec- 
centricities and induces secular resonances on the icy gi- 
ants. There is a period of crossing orbits which eventually 
scatters the icy giants outwards, which enter the planetesi- 
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mal disc (|Levinson fc Morbidelli|[2003l ). The icy giants eject 
inwards some planetesimals to Jupiter and Saturn, which 
makes them migrate even further away. The planets finish 
their migration once the planetesimal disc is completely de- 
pleted. The final configuration of the first Nice model is close 
to the observed values for both the semi-major axes an d 
the mean eccentricities of the pla nets l|Tsiganis et al.ll2005l ). 
Recently. IMorbidelli etafl i|2009l ) extended their original re- 
sults showing that encounters of Saturn with an ice giant 
lead to the correct secular evolution for the eccentricities of 
Jupiter and Saturn, and not the passage through the 1:2 
mean-motion resonance, as it was originally proposed. 

An important open issue in the Nice model is related 
to explain the initial conditions assumed in the model, 
from earlier stages of the formation. This question has been 
recently investigated, assuming an earlier multi-resonant 
configuration for the already formed giant planets, be- 
sides pushing the planets into such configuration, by in- 
trod ucing additionally semi- major axis and eccentricity de- 
cay l|Batygin fc Brownll2010l ). Yet, the issue remains to ex- 
plain those initial conditions. 

The aim of this paper is to study the statistics of the 
final configurations of planetary systems for some models 
of planetary accretion and evolu tion, which ar e based on 
a model originally introduced by lLaska r (2000a). The final 
orbital parameters we obtain could be used as the initial 
conditions for detailed calculations on the planetary evolu- 
tion, e.g., aimed to describe the architecture of the planetary 
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systems formed in the spirit of the Nice model, or to address 
the dynamical stability of these systems in statistical terms. 
Yet, it is not our purpose to obtain an accurate comparison 
with the initial conditions of the Nice model, or with the 
architecture of the Solar System. Our aim is to understand 
the consequences of incorporating in a simple form some 
important physical effects on the accretion processes. The 
variations of Laskar's model included account for the sys- 
tematic conservation of the total angular momentum of the 
system, and incorporate mass- and velocity-selection mecha- 
nisms that tune the accretion probability; note that Laskar's 
model actually conserves the angular momentum, but not 
the methodological procedure used to achieve a fixed final 
angular momentum deficit (AMD). These aspects were not 
considered in Laskar's model, though they constitute impor- 
tant building blocks of the current planetary formation the- 
ories l|Safronovlll972l ; |Polack et al.lll996l ; lArmitage|[201ol '). By 
comparing the outcome of the mechanisms that we include, 
we are able to relate them, in a statistical sense, to certain 
effects manifested by the final configurations. This is done 
directly by correlating the masses and semi-major axes of the 
planets formed in a large number of realizations of the mod- 
els. Our results indicate a power-law behaviour of the mass 
in terms of the semi-major axis, whose exponent depends on 
the different models, and on the initial linear mass-density 
distribution assumed. In particular, we show that includ- 
ing both mass- and velocity-selection mechanisms leads to 
the systematic appearance of planets in certain specific lo- 
cations. Moreover, by introducing properly-scaled variables, 
we show that these results are universal with respect to the 
total angular momentum of the system, the mass of the plan- 
etesimal disc and the mass of the central star. 

The paper is organised as follows: In Section [2] we de- 
scribe Laskar's simplified model of planetary evolution and 
accretion in detail, clarify its predictions, and outline a short 
critique aimed to motivate the variations that we incorpo- 
rate in the model. Specific details of the variations we imple- 
ment and the corresponding numerical results are described 
in Section [3] In Section [4] we address how variations of the 
total angular momentum of the system, the mass of the plan- 
etesimal disc, and the mass of the central star affect our 
results. Finally, in Section [5] we summarise our results and 
conclusions. 
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Here, Q is the gravitational constant, n is the mutual dis- 
tance between the i-ih planetesimal and the central star, 
Tij is the distance between the i-th and j-th planetesimals, 
Po = — 2 =i an d Pi are the canonical momenta as- 
sociated with the star and the planetesimals, respectively. 
Clearly, Ho corresponds to N independent two-body Kepler 
problems among each planetesimal and the star, and Hi 
is the perturbing term, which includes the mutual gravita- 
tional interaction among the planetesimals and the indirect 
term. Laskar's model is formulated for the spatial problem; 
the numerical simulations are illustrated for the planar case. 
From now on, we focus on the latter case considering copla- 
nar orbits. 

The dynamics of Ho is completely integrable: For each 
planetesimal the energy Ei — — /im»/(2aj) and the barycen- 
tric angular momentum U = mi\pa,i(l — e 2 )] 1 ^ 2 (/x = QMo) 
are conserved quantities. This implies the conservation (for 
Ho) of the semi-major axis a< and eccentricity e, of the ellip- 
tic orbits of the planetesimals. The non-integrable character 
of H is manifested in the long-term evolution and is due to 
the per turbing terms of Hi. To model these non-integrable 
effects, iLaskarl l|2000ah considers the secular system, where 
the equations of motion are averaged over the mean longi- 
tudes. This simplification is tantamount to excluding effects 
related with mean-motion resonances. 

Firstly, considering the case without collisions, the av- 
eraged system conserves the energy of each planetesimal, 
but exhibits a slow chaotic diffusion of the individual ec- 
centricities ei. This chaotic diffusion is constrained by the 
conservation of the total angular momentum of the system 
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Including the conservation of the energy, this leads to the 
conservation of the total angular momentum deficit (AMD) 
of the system. The latter is expressed as 
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2 A SIMPLIFIED MODEL OF PLANETARY 
FORMATION AND EVOLUTION 

2.1 The model 

Lask ar's simplified model of planetary formation and evolu- 
tion l|Laskanl2000al lbh considers a system which consists of 
one large central body of mass Mo and a large number of par- 
ticles of small mass rrii, i = 1 . . . N (N 3> 1), that interact 
under their mutual gravitational attraction. The Hamilto- 
nian of this (N + l)-body problem, in a heliocentric frame, 
can be written as H = Ho + Hi , where Hq and Hi are given 
by 
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Physically, the total AMD of a planetary system is a mea- 
sure of the circularity of the planetary orbits, since for small 
eccentricities C; ~ niilpUH] 1 Conversely, large values of 
Ctot indicate the possibility of planetary collisions. 

Pair- wise collisions in this simplified model are assumed 
to be totally inelastic and lead to accretion. Assuming that 
mass and linear momentum are conserved in the collision, 
rrii^j — rrii + rrij and P;ej = Pi + Pj , the conservation of the 
angular momentum follows. From these conservation laws, 
the orbital elements of the accreted particle i © j can be 
calculated: First, the semi-major axis of the accreted parti- 
cle is obtained from its energy, which can be expressed as 
Ei®j = Ei + Ej + SEij. The change in the energy is given 
by 
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Then, using the conservation of angular momentum of the 
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colliding particles yields the eccentricity of the accreted par- 
ticle. 

As a crucial final element, lLaskarl (|2000bh demonstrates 
that under this accretion scheme, the local AMD of the ac- 
creted particles satisfies 



(6) 



This fact motivate s the so-called AMD stability crite- 
rion (Laskar 2000a): A planetary system is AMD stable if 
its total AMD is not sufficient to permit planetary collisions. 
This criterion suffices to ensure the long-time stability of the 
averaged system, since new collisions are not possible; how- 
ever, this is not necessarily true for the complete system 
H , since the latter includes effects linked with short-period 
resonances which could still induce collisions. 

Algorithmically , the model is implemented as fol- 
lows l|Laskarll2000ar ). Initially, the semi-major axes are dis- 
tributed homogeneously throughout the disc and their as- 
sociated masses are fixed by the linear mass distribution 
p a (a); likewise, the initial eccentricities are drawn from 
the distribution p s (e). The planetesimals are labelled in 
non-decreasing order according to their semi-major axis: 
cii < a 2 < ... < ajv, where is the initial number of 
planetesimals. Note that no angular variables are specified 
so far. The iterations proceed as follows: One planetesimal i 
is chosen at random from the list of planetesimals together 
with one of its neighbours, say i + 1. If the co-focal ellipses 
associated with their trajectories display an intersection -in 
the purely geometrical sense-, which is expressed as 
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then the particles are accreted at one of the intersection 
points. Note that the relative orientation of the orbits is 
a random variable defined in an appropriate interval, which 
guarantees the existence of a geometrical intersection; the 
necessary condition Eq.[7]thus corresponds to anti-alignment 
(8 = 7r) of the ellipses. The orbital elements of the accreted 
particle are calculated and the list of planetesimals is up- 
dated (as mentioned above, the AMD of the accreted par- 
ticle is smaller than the sum of the AMDs of the colliding 
particles). Between collisions, a random walk in the space 
of eccentricities is implemented, which models t he ch aotic 
secular evolution of the system ( Laskarl 1 19891 . I1990T ). re- 
stricted by the conservation of the total angular momentum. 
These steps are iterated until the total AMD of the system 
is smaller than a critical value C cr , thus ensuring that no 
planetary collisions are possible. Clearly, the final planetary 
system formed is the result of two competing processes, a 
chaotic diffusion due to the secular evolution of the averaged 
system that promotes collisions, and the circularisation due 
to accretion that inhibits them. 

An important simplification of Laskar's model relies on 
the introduction of independent stochastic discrete time- 
steps or iterations, instead of computing the detailed dynam- 
ical evolution of the many-body system, which is a compli- 
cate and time-consuming task. The consequence of this sim- 
plification is the limitation to provide any physically relevant 
time-scale. Indeed, each iterate of the algorithm described 
above involves the secular chaotic diffusion in the space of 
angular momentum, and the computation of the resulting 
accreted particle from a total inelastic collision of two plan- 
etesimals. All details of the actual physical processes are 



abandoned by modelling them as a Brownian motion in the 
space of eccentricities (or equivalently in the space of angu- 
lar momentum), and by taking a random orientation of the 
elliptical orbits and one of its geometrical intersections. Ad- 
ditionally, each time-step is computed independently from 
the previous one, in the sense that there is no memory. 
Each of the physical processes involved in each iteration have 
different time-scales associated, which depend on the semi- 
major axis, the local mass-density, etc. This implies that 
two distinct discrete time-steps may involve quite different 
scales of the physical time. Therefore, an important limita- 
tion of Laskar's simplified model -which is inherited by our 
variations-, is the impossibility to answer questions on the 
physically relevant time-scale based on the simulations. 

Despite of the lack of a characteristic time-scale, the 
simplicity of Laskar's model is reflected in the fact that 
it permits certain analytical treatment that yields concrete 
predictions. These results follow from the minimum possi- 
ble AMD value of two colliding particles, that takes place 
with the outermost located at perihelion and the in nermost 
at aph elion. After some algebraic manipulations (|Laskad 
|2000bJ), the minimum or critical AMD C cr , which allows 
for that collision, can be obtained. Using the asymptotic 
properties of C cr properly scaled, Laskar obtains scaling 
laws fo r the spacing between adjacent planets and their 
masses (Laskar 2000a). We write Laskar's results as 
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Here, S a = ai+i—ai denotes the spacing between two nearby 
planets (Sn — 1), p a {a) is the initial mass distribution of 
planetesimals, m(a) is the mass of the planet, and k is a 
dimensionless constant which depends on the ratio of the 
colliding masses. Equation [8] provides the spacing of the 
planets, while Eq. [9] provides their masses; both are given in 
terms of the initial mass distribution p a {a). We notice that 
Ccr in these equa tions is actually an unknown quantity. 

lLas kar ( 2000a|) illustrated his findings numerically using 
a constant linear mass-density distribution p a (a) = Co f° r 
the planetesimals, a total mass of the disc rtlT = 8 x 10~ 6 
(in solar masses) and for the AMD the value Cfi Ila i = 16 X 
10 -8 as the statistically significant quantity. His numerical 
results, in terms of the number of the planet n, showed good 
agreement with respect to the analytical predictions, given 
by Eqs.HandU 



2.2 Comments 

Laskar's model described above is elegant for its simplic- 
ity and predictions, regardless of the lack of any char- 
acteristic physical time-scale. In deed, Eq. [8] c an be inte- 
grated explicitly in many cases dLaskarl |2000al ). In partic- 
ular, for p a (a) = (, a CL Za , we have that for z a 7^ —3/2 



za/3+1/2 



this yields a power-law spacing distribution a 
a z a /3+i/2 k n ^ an example, z a = yields the square- 
root behaviour al/ 2 ~ n. Likewise, for z a = —3/2 it yields 
a n = ap exp \Ccr 3 P~ 1 ^ 6 n/(kC^ 3 )] , which is basically the 
usual Titius-Bode law (Nieto 1972). Similar expressions can 
also be obtained for the mass in terms of n using Eq. [5] 

While the results are attractive, there is an important 
subtle point in Eqs. [8] and [9] as they stand. As mentioned 
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Figure 1. Frequency distribution of Cfl na i/Ltot for 4096 inde- 
pendent realizations of Laskar's model keeping fixed the total 
angular momentum of the system. In these calculations we fixed 
Ltot = 2.2212 X 1O~ 2 M0 AU 2 yr _1 , which is roughly the angular 
momentum of the planets of the Solar System. 



there is no physical constraint with respect to the relative 
velocity o f the coll i ding particles, as stated by Safronov's 
criterion (|Safronov|[l97a ). The latter establishes that ac- 
cretion is possible if the motion of the colliding particles 
satisfies v Ic i ^ « csc , where v re i is the relative velocity of 
the colliding particles at the point of collision, and v csc is 
the escape velocity of the local two-body problem. The for- 
mer depends -among other parameters defining the elliptic 
orbits- on the relative orientation 8, while the latter can be 
written as f 2 sc = 2/i/rn, where m is the Hill radius de- 
fined with respect to the dominating mass. Moreover, the 
accretion probability in Laskar's original model is indepen- 
dent of the orbital and physical parameters of the parti- 
cles, bei ng in this sense tantamou nt of an orderly type of 
growth (|Wetherill fc Stewart! [1989). While this is certainly 
a good starting point, it is intuitively important to include 
a mass-dependence in the accretion probability, since it af- 
fects the collision cross-section. A particular case of inter - 
est is the so-called runaway growth l|Greenberg et al"] [l978'). 
during which the more massive bodies grow faster than the 
lighter ones. 



above, the total AMD is constant between consecutive col- 
lisions and decreases only when they occur. Then, below a 
critical value C cr it remains constant. Note that C cr is not 
a constant of motion in the dynamical sense. For the com- 
parison with the numerical simulations, Laskar uses the sta- 
tistically significant quantity Cfinai, the final AMD. While 
Cflnai is not a constant of motion, because different statisti- 
cal realizations of the model yield different values, the nu- 
merical results are quite satisfactory. Indeed, in order to 
achi eve the fixed v alue Cfinai, as a methodological proce- 
dure lLaskarl (|2000al ) excites the eccentricities of the planets 
at the end of the simulation, thus losing the constancy of 
the total angular momentum of the system; note that it is 
here that Laskar's implementation fails to conserve the to- 
tal angular momentum of the system. This observation is 
important since the conservation of angular momentum is a 
corner stone of the theories o n planetary formation and evo- 
lution (|Fernandez fc Irj||l984T ). Therefore, for the integration 
of Eqs. [8] and [9] some previous knowledge of the distribu- 
tion of the Cfinai is required, or at least its dependence with 
respect to other initial parameters, e.g., p a (a). Note that 
knowledge of the distribution of Cfi na i does provide means 
to consistently define C cr . 

To illustrate this point, in Fig. Q] we plot the frequency 
distribution of the final AMD normalised to the total angu- 
lar momentum of the system, Cfi na i/Ltot, for a large number 
of realizations of Laskar's simplified model, keeping the to- 
tal angular momentum of the system constant. The figure 
clearly illustrates the lack of constancy of the final AMD 
for different realizations. We further observe that C cr can 
be defined as the supremum of Cfinai; while this is a con- 
sistent definition, the numerical calculation of such quantity 
may be involved due to the poor statistics at the tail of the 
distribution. 

Aside from the lack of conservation of L to t , in the nu- 
merical implementation of Laskar's model, all relative ori- 
entations of the elliptic trajectories of the colliding particles 
that intersect lead to accretion with equal probability. In 
contrast, in the derivation of Eqs. [8] and [9] Laskar imposes 
a very specific collision among the particles. In either case, 



3 VARIATIONS ON LASKAR'S MODEL 

In this section we describe the variations and the numer- 
ical results obtained by implementing some variations to 
Laskar's simplified planetary model. These variations im- 
pose, firstly, the conservation of the total angular momen- 
tum of the system, and secondly, address the consequences 
of using mass-dependent accretion processes and velocity 
selecti on mechanisms, not considered previously (|Laskarl 
without leaving the overall simplicity in the for- 
mulation and implementation of the model. 

3.1 General characteristics of the simulations 

For the numerical results presented in this section, we con- 
sider that the mass of the host star is one solar mass 
(Mo = IMq), which we shall use as the unit of mass. We 
shall also fix the astronomical unit (AU) as the unit of dis- 
tance and the yr as the unit of time. In these units we have 
H = GM = 4tt 2 AU 3 yr _2 MQ 1 . In addition, we shall fix 
the total mass of the disc to M disc = 1.3413 x 10" 3 M Q 
and the total angular momentum of the system to Ltot = 
2.2212 x 10~ 2 M© AU 2 yr _1 . These values are roughly the 
current values of mass and total angular momentum of the 
planets in the Solar System, respectively. These parameters 
define the physical quantities of our simulations. 

In order to define the initial conditions, we must spec- 
ify the form of the initial linear mass-density distribution 
Pa{a) (or equivalently the initial surface density p a (a)/a) 
and the initial eccentricity distribution p e (e). These distri- 
butions and the extension of the disc are constrained by the 
total mass and the total angular momentum of the system. 
In the continuum limit, the total angular momentum of the 
system is written as 

/* a max /* e max 

Ltot = / / pa(a)p e (e)y / pa(l - e 2 ) dade. (10) 

Assuming that a m i n and e max are given quantities, Eq. [10] 
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can be used to define the maximal extension of the ini- 
tial disc, ttmax • Notice that Eq. \W\ can be explicitly inte- 
grated for a homogeneous eccentricity density p e (e) = £ e 
and a power-law linear mass-density distribution, p a (a) = 
( a a Za |Hernandez-Menallunpublishedh . 



F or a fair comparison with Laskar's results (|Laskarl 
2000a), our simulations are initiated with a large num- 
ber (typically 10000 bodies) of equal mass planetesimals, 
which is the most constraining case in terms of equivalence 
of th e planetesimal location jNamouni. Luciani fc Pellatl 
Il996h ; similar results are obtained for a larger initial num- 
ber of planetesimals. We also considered the planetesimals 
on coplanar orbits, defining their initial semi-major axes and 
eccentricities from a homogenous linear mass-density distri- 
bution p a (a) = Co an d a homogeneous eccentricity distri- 
bution p e (e) — Ce, respectively. Using a m j n = 0.1 AU and 
e m ax = 0.3 in Eq. [I0]yields for the maximal extension of the 
disc a max = 15.93 AU. This value of a max is smaller (by a 
factor 2) than the 30 — 35 AU estimate for the outer edge of 
th e original planetesimal disc of the Solar System, proposed 
bv lLevinson fc M orbidclli (2003]). However, this value is fully 
consistent with the range for the initi al conditions used fo r 
the giant planets in the simulations of lTsiganis et al.l l|2005h . 
The contributions of the massive disc (30 — 50 M§) that 
these authors include beyond the orbits of the planets up 
to 30 — 35 AU are excluded from the values Mdi sc and Ltot • 
Notice t hat the outer disc ha s also a constant linear mass- 
density (|Tsiganis et al.l [2005 V Finally, we emphasize that 
the use of Eq. [10] to define the maximal extension of the 
initial disc, a max , is yet another difference with respect to 
Laskar's implementation. 

It is also worth describing the implementation of the 
Brownian motion that accounts for the chaotic evolution 
of the averaged equations. Instead of implementing it in the 
space of eccentricities, we accomplished this by the exchange 
of angular momentum between a pair of arbitrary planetes- 
imals. Denoting the initial angular momenta of the chosen 
planetesimals by h and l 2 and the corresponding final ones 
as I'x and 1' 2 , we write l[ — h + SI and 1' 2 = l 2 — SI. This 
ensures the conservation of the total angular momentum. 
The angular momentum exchanged is written as SI — 
Here, /3 ^ 1 is a positive dimensionless parameter related 
to the diffusion constant; in our simulations we considered 
the value /? = 10~ 3 which i s small enough to model the 
slow secular chaotic diffusion (Laskar 1989, 1990). The ran- 
dom walk character is provided by the stochastic variable £, 
which is a uniformly distributed random number in the in- 
terval [— c' 2 , Ci]. Here, c' x and c' 2 axe the largest positive num- 
bers (with units of angular momentum) that simultaneously 
satisfy the inequalities c[ ^ C\, c[ ^ l 2 , c' 2 ^ C 2 , c' 2 ^ 1%, 
where C\ and C 2 are the AMD's of the planetesimals. This 
definition of 51 ensures that both l[ and 1' 2 lie in the cor- 
rect intervals and the new eccentricities are well-defined. 
Note that this definition permits collisions of very eccentric 
(d ~ 1) planetesimals with the host star |Weidenschillinel 
1 19751 ). This is consistently taken into account in our simu- 
lations whenever the perihelion of a particle is smaller than 
the radius of the Sun (4.65 x 10" 3 AU). 



3.2 Laskar's model including angular momentum 
conservation: orderly growth 

We shall discuss firstly the case where we impose con- 
sistently the conservation of the total angular momentum 
within the simplified model of Laskar. This is aimed to clar- 
ify the differences that arise with the variations we imple- 
ment with respect to Laskar's results. As mentioned above, 
the simulations end once the system satisfies the AMD sta- 
bility criterion. We carry out accretion and orbit evolu- 
tion by selecting at random a planetesimal and one of its 
neighbours for accretion, and then proceed with the ex- 
change of angular momentum between (many) pairs of plan- 
etesimals, also chosen at random. The mass distribution 
evolves smoothl y, which makes this case similar to an or- 
derly growth (cf. Armitagc 2010). 

In Fig. [2] we present the final planetary configurations 
by combining the results of all 4096 realizations of the 
model, plotting the mass and eccentricity of the formed 
planets in terms of their semi-major axis. Figure [3] illus- 
trates the resulting final configuration of one planetary 
system taken at random from the set of realizations. In 
both figures (and throughout this section), the final plan- 
etary masses are given in units of the mass of Jupiter Mj 
(1M Q = 1047.56Mj; lMj = 317.83M®). 

The results presented in Fig. [2ja) show an increasing 
trend of the planetary masses up to a « 14.5 AU, where 
a rapid decrease of m in terms of a takes place. The lat- 
ter is due to the finite size of the planetesimal disc and the 
conservation of the angular momentum, and corresponds to 
the location of the outermost planet of the planetary sys- 
tems formed. We observe that the largest planetary mass 
attained is ~ 0.2Mj, which is a comparatively small value 
with respect to the mass of Jupiter. This is a consequence 
of the large number of planets formed in each simulation, 
producing an average of 30.9 ± 1.7 planets (the error is the 
standard deviation). 

As illustrated in Fig. individual planetary systems 
may not display the purely initial increasing behaviour of 
m(a) observed in the data of the ensemble. The results of 
the ensemble are statistical and it is thus m eaningful to de - 
fine the average behaviour of the ensemble. Laskar (2000a) 
computes the average of the final planetary mass and semi- 
major axis for fixed planet number n. Since n is not an 
observable quantity, we have opted to average the results 
over a small interval of the semi- major axis, which was fixed 
to 0.2 AU. The resulting (binning) average is represented by 
the star symbols in Fig. [2] For comparison with Eq. [9l we 
fit the averaged data with the power-law function 

m(a) = Aa". (11) 

According to Eq. [9] we should expect v = 1/2. The contin- 
uous curve in Fig. [2Ja) illustrates the resulting best-fitting, 
which in this case was calculated up to a — 14.5 AU where 
the decreasing branch of the last planet sets in; the best- 
fitting yields v « 0.59. This value of v is close but larger 
than Laskar's predicted value 1/2. We attribute this differ- 
ence to the (implicit) dependence of the distribution of C cr 
on the specific parameters of the model; this dependence is 
not considered in the integration of Eqs.[8]and[9] In Fig.^a) 
we also plot the best-fitting curve, in order to compare the 
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Figure 2. (a) Final planetary mass and (b) eccentricity in terms 
of the corresponding semi-major axis for 4096 realizations of 
Laskar's simplified model, imposing the conservation of the total 
angular momentum of the system. Each dot represents a planet 
formed in one of the simulations. The star symbols represent local 
averages taken over a semi-major axis interval of 0.2 AU. In the 
top panel, the curve represents the power-law fit (Eq. 1116 of the 
averaged data in the interval a £ [0, 14.5]. 
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Figure 3. (a) Final planetary mass and (b) eccentricity in terms 
of the semi-major axis for a typical (randomly chosen) plane- 
tary system, obtained by implementing Laskar's simplified model 
with conservation of the total angular momentum (orderly-growth 
model). The continuous curve in the top panel is the the power- 
law fit (Eq. 1 1 1 D obtained from the averaged data of the ensemble 
of Fig.rja). 



results of a specific random planetary system with the re- 
sults averaged over the ensemble. 

Figure [2jb) displays the results of the final planetary 
eccentricity for the ensemble on a logarithmic scale, and 
Fig. [3{b) the corresponding results for a typical member of 
the ensemble. We observe a discontinuous behaviour around 
a ~ 0.1 AU. This feature is similar to the discontinuity re- 
lated to the location of the outermost planet discussed in 
the mass diagram; in this case, it indicates the localisation 
of the innermost planets. The figure shows that th e final 
eccen tricities are smaller for larger semi-major axis l|jonesl 
I2004T ). Planets with smaller masses display larger eccentric- 
ities in comparison with planets with larger masses. This is 
a consequence of the accretion processes, which tend to cir- 
cularise the planetary orbits, consequently causing the more 
massive planets to exhibit smaller eccentricities than those 
less massive. 



3.3 Mass-selection mechanism: runaway growth 

As mentioned earlier, in Laskar's original model the accre- 
tion probability is independent of the mass of the planetes- 
imals. Yet, mass increases the collisional cross-section and 
thus promotes accretion, at least once the gravitational at- 
traction among planetesimals starts to dominate their dy- 
namics. We consider now the inclusion of such effects in 
the model. We address the case where the massive bodies 
grow quite fast due to their mass, and thus increase their 
difference in mass with respect to the lighter bodies. This 
happens until they exhaust the available mass in their sur- 
roundings. In what follows, we shal l refer to this process 
as mass-selection or runaw ay growth (jGreenberg et al . 1978; 
IWetherill fc StewartJl 1989). Notice that, while we are includ- 
ing now the mass to select the particle that accretes, the rel- 
ative velocity of the colliding particles is not yet taken into 
account, which is anoth er important pro perty that enhances 
gravitational focussing (|Armit agc 2010). 

We implement the mass-selection growth as follows. At 
each iteration of the model we define a threshold mass m cut , 
which is a uniformly-distributed random number between 
zero and the mass of the heaviest planetesimal. We then se- 
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lect at random a planetesimal whose mass is larger than this 
threshold mass, and implement the accretion process with 
one of its randomly-selected neighbours; the relative orien- 
tation of the elliptic orbits of the colliding particles is con- 
sidered as before. With this simple implementation, once a 
particle has a mass slightly larger than the rest, i.e., from the 
first iteration, the probability that this particle is selected 
for collision in the next iteration is dramatically increased 
in comparison to other particles. This promotes accretion 
of the more massive particles, and hence mass growth in a 
runaway-like form. Additionally, it will systematically open 
a gap around the location of the accreted particle, since the 
distance to the nearest neighbours increases. Eventually, the 
nearest particles will be far enough to avoid collisions with 
this particle, i.e., the geometric condition [7] is not satisfied. 
In this case the planetesimal becomes locally isolated, un- 
til another particle is close enough again or enough angular 
momentum has been exchanged to allow for a new collision. 
If the particles that satisfy the condition rrii > m cut can not 
collide with their neighbours (the isolation mass is locally 
reached), we then reset the value of m cu t to a smaller value, 
and proceed as described above. Clearly, this implementa- 
tion promotes accretion of the more massive planetesimals. 
The chaotic diffusion following the accretion processes is im- 
plemented as before. Again, the simulations are iterated un- 
til the condition of AMD stability is reached. 

In Fig. |4] we present the results of the runaway (mass- 
selection) growth model, illustrating the combined results of 
the masses and eccentricities of the final planetary configura- 
tions in terms of their semi-major axes, for 4096 independent 
realizations of the model. Figure[5]displays similar results for 
a single planetary system chosen at random. As illustrated 
by the results of the ensemble, the mass increases with the 
semi-major axis up to a ~ 13.4 AU, where a sudden decrease 
of the mass takes place. This is similar to the case of orderly 
growth, where the decreasing branch marks the location of 
the outermost planet formed. As done before, we fit the local 
averaged mass (up to 11.4 AU) with the power-law function 
(Eg. Illj l. This yields the exponent v « 0.66, which is larger 
than the one obtained for the averaged data of Fig. [2ja). 
We have chosen to fit up to a = 11.4 AU because, beyond 
this value, the averaged mass displays a strong oscillation. 
Notice that the overall mass-scale is enhanced by a factor 
close to 2 with respect to the orderly-growth case, but the 
largest accreted particles are at most 0.4Mj. The average 
final number of formed planets is reduced to 14.6 ± 1.15. 

It is interesting to compare the density of points in the 
mass-semi-major-axis diagram. For a given small interval of 
a, in Fig. [2ja) the density of the ensemble of planets is quite 
uniform with respect to the mass they reach (except towards 
the upper edge). This statement holds essentially for all a be- 
fore the last-planet branch appears. On the other hand, the 
density of points in Fig. |4ja) is manifestly non-uniform, at 
least for semi-major axis larger than ~ 2 AU. Indeed, consid- 
ering a small interval of semi-major axis, there is a marked 
clustering of points towards the upper and lower edges of 
the mass. According to these observations, this runaway- 
growth mechanism yields larger final planetary masses, and 
thus a smaller number of planets, but also a systematic oc- 
currence of small-mass planets for essentially all values of a. 
Since these planetary systems are AMD stable, these small- 
mass planets are far from planetary collisions. In this case, 
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Figure 4. Same as Fig. f2] combining 4096 realizations of 
the model, when the implementation includes runaway (mass- 
selection) growth. In the top panel, the fit with the power-law 
Eq. llll of the averaged data, represented by the continuous curve, 
was carried in the interval a S [0, 11.4]. 



possible mean-motion resonances not included in the model 
could increase the ir eccentricities and e ventu ally yield new 
collisions (see e.g. lLaskar fc G astincau (2009])). 

The eccentricities of the planetary ensemble obtained 
in this case, illustrated in Fig. |4jb), differ also from those 
resulting from orderly growth, cf. Fig. [2jb). Figure |4jb) 
shows the lack of uniformity in the density of points; yet, 
we note there is a density concentration towards compar- 
atively smaller values of the eccentricity. That is, compar- 
ing the average value of the eccentricity for a given interval 
around a for orderly and runaway growth, in the latter the 
average eccentricity is smaller, except for large values of a. 
Therefore, the mass-selection mechanism induces more cir- 
cularised orbits of the planets, and in this sense, the final 
configurations are more stable. We also note that in the 
present case, individual planets typically with small a, may 
reach comparatively larger values of the eccentricity. 



3.4 Velocity-selection mechanism 

Another important factor that influences the accretion prob- 
ability is the ratio between the relative velocity of the col- 
liding bodies and their escape velocity. In particular, this 
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Figure 5. Same as Fig. [4] for one typical planetary system ob- 
tained in the runaway (mass-selection) growth model. The con- 
tinuous curve in the top panel corresponds to the power-law fit 
of the averaged data in Fig. |4ja) . 
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Figure 6. Same as Fig. [2]for the combined results of 4096 simu- 
lations of the velocity-selection accretion model. In the top panel 
the continuos curve represents the fit of the power-law carried out 
in the interval a £ [0, 9.34]. 



quantity enters in the collisi onal cross-section through the 
gravitational focusing factor (|Wetherill fc Stewart||l989T ) 

Fg = l+t&c/w?ei- (12) 

Here, v rc i is the relative velocity of the colliding particles 
(at infinity), and v esc is the escape velocity in the local two- 
body problem. For w re i/«esc <C 1 the gravitational factor 
enhances the collision cross-section. Moreover, in this case 
the colliding bodies are expected to be bound toget her and 
thus promote planetary growth fsee lArm itagc 2010l). 

One naive way of implementing these effects into 
Laskar's simplified model with the conservation of angu- 
lar momentum, is by selecting the relative orientation of 
the ellipses of the colliding bodies, such that the condition 
v IC \/v csc <C 1 is satisfied. That is, we allow accretion of par- 
ticles only when Safronov's condition is fulfilled. Carrying 
on these simulations resulted in systems where the AMD 
stability condition was not fulfilled, in general. Therefore, 
the elliptic trajectories of nearby bodies may intersect and 
thus collide, but for those collisions no value of 9, the rel- 
ative orientation of the ellipses, exists that v Te \/v csc <C 1 is 
fulfilled. 

We therefore relaxed Safronov's condition for accre- 
tion, and fixed 6 as the relative orientation of the ellipses 
where v Te \ achieves its minimum value. By doing this, we 



promote a moderate enhancement of the collisional cross- 
sections by the gravitational focussing factor, as well as 
quasi-tangent ial collisions. T he latter is physically plausible 
for accretion ((Ohtsuki 199^). We note that with this imple- 
mentation, the change in the absolute value of the energy 
after accretion is minimal, cf. Eq. [5] thus minimising the 
net migration of the accreted particle. We shall refer to this 
mechanism as the velocity-selection mechanism. 

The results of these simulations are shown in Fig. [6l for 
the combined results of the ensemble, and in Fig.[7]for an in- 
dividual planetary system taken at random. As shown in the 
results for the ensemble, the mass increases with the semi- 
major axis up to a Ri 12.6 AU, where the feature associated 
with the location of the outermost planet appears. Fitting 
Eq. [TT]with the local average mass up to 9.34 AU yields the 
power-law exponent v m 0.74. This value of v is even larger 
than the one obtained for runaway accretion. The mass-scale 
is larger as well, being enhanced by about 25% with re- 
spect to the mass-selection mechanism; the maximum mass 
reaches approximately 0.53 Mj. Correspondingly, the aver- 
age final number of planets is further reduced to 9.3 ± 0.7. 
As done before, the limit for fitting the power-law was con- 
sidered due to the appearance of the strong oscillations in 
the averaged data, as illustrated in Fig. |6ja). 
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Figure 7. Same as Fig. [6] for a random planetary system ob- 
tained in the velocity-selection model. The continuous curve in 
the top panel is the power-law fit of the averaged data displayed 
in Fig. Oa). 
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Figure 8. Same as Fig. [2]for the combined results of 4096 simu- 
lations of the oligarchic growth model. In the top panel the con- 
tinuos curve represents the fit of the power-law carried out in the 
interval a £ [0, 8.6]. 



The density profile displayed in Fig. [6ja) resembles the 
one obtained for the orderly growth model, being more uni- 
form than the density profile for runaway growth. Yet, there 
is a qualitative difference with respect to the former results, 
which is manifested by a gap in the mass-semi-major axis 
diagram. Indeed, there is a region in the m vs a diagram, 
defined beyond 10.7 AU (for small masses), where no planet 
is formed. That is, planets formed beyond ~ 10 AU have a 
mass above a certain lower bound, which depends on a. Ac- 
tually, a similar feature can also be observed for small values 
of a. We also note the less-massive planets formed at inter- 
mediate semi-major axis are consistently more massive than 
the corresponding ones for the mass-selection mechanism. 

With respect to the eccentricities of the ensemble, as 
illustrated in Fig. [IJb), we confirm the overall resemblance 
with the results obtained for orderly growth, where 9 is a 
uniformly-distributed random variable in the proper inter- 
val. In the present case, the distribution of final planetary 
eccentricities yields the possibility of larger values of the ec- 
centricity. This can be noticed, e.g., in the branch of small 
a associated with the location of the innermost planets, and 
for large values in the sudden increase of the average eccen- 
tricity of the outermost planets. 



3.5 Mass- and velocity-selection mechanism: 
oligarchic growth 

We consider now the combined effect of the mass- and 
velocity-selection mechanisms. To this end, we select a par- 
ticle according to the mass-selection mechanism described 
previously, and then select the relative orientation 9 of the 
elliptic orbits of the colliding particles such that v IC \ is a 
minimum. This case is somewhat a nalogous to oligarchic 
growth models (|Kokubo fc IdalFl99ct ). where after an initial 
runaway phase, gravitational focussing is enhanced by all 
the dominating protoplanets but on a slower time-scale. We 
shall thus use this name to refer to the present case. The 
combined results of 4096 simulations are shown in Fig. [8] 
and the result of one randomly selected realisation in Fig. [5] 

As illustrated in Fig.[8ja), the results display an increas- 
ing behaviour of the mass upon the semi-major axis, until 
the feature associated with the position of the outermost 
planet appears. In this case, fitting the averaged data with 
Eq. [TT]in the interval a € [0, 8.6] yields v ^ 0.58. This value 
is comparable (marginally smaller) than the value obtained 
for orderly growth. Despite this, the overall mass-scale is 
slightly larger than in previous cases, reaching the value 
0.63 Mj, consequently yielding a smaller average number 
of formed planets: 6.9 ± 0.7. 

With regards to the final eccentricities, Fig. |8jb) illus- 
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Figure 9. Same as Fig. [8] for a randomly selected planetary sys- 
tem obtained in the oligarchic growth model. The continuous 
curve in the top panel is the power-law fit of the averaged data 
displayed in Fig.[§Ia). 



trates the enhancement of the final eccentricity of the orbits. 
Indeed, we observe an increased variance of the final eccen- 
tricity at all scales of the semi-major axis, as well as their 
average value. In particular, we note that for the outermost 
planet, the corresponding eccentricity may reach values close 
to 0.2. Interestingly, this indicates that the distribution of 
Cflnai has larger values than in the other cases. 

The density profiles displayed in the mass and eccentric- 
ity diagrams are interesting and differ clearly in the struc- 
ture from the other cases. Indeed, not only the location of 
the outermost and innermost planets is noticeable, but actu- 
ally the location of other planets can be estimated from the 
darken clumps that appear in the figures. This can also be 
observed in the oscillations displayed by the averaged data, 
which, in comparison with the cases discussed before, are 
more prominent and appear even for smaller values of a. 

We illustrate these observations in Fig. [TUJ In Fig. HOf a) 
we plot the number of final planets h(a) (on a log scale) lo- 
cated in an interval of size Aa ~ 0.2 AU centred around 
a semi-major axis a, for the four different growth models 
discussed. The dashed curves are the best-fitting power-law 
corresponding to each data set, which we denote hat{a). In 
all growth models, we observe a clear excess of particles 
(with respect to the best-fitting) towards the outer edge of 
the planetary systems, which is associated with the location 




Figure 10. (a) Number of planets h(a) formed in an interval of 
size Aa ~ 0.2 AU around a. The symbols identify the different 
growth models: Triangles (red) for orderly growth (Laskar's model 
with conservation of Ltot), squares (cyan) for runaway growth, 
empty circles (magenta) for the velocity-selection mechanism, and 
filled circles (blue) for oligarchic growth. The dashed lines are the 
best-fitting power-law hfi t (a) for each case, (b) Relative residuals 
with respect to /ifl t (a). The peaks indicate a systematic enhance- 
ment of the formation of planets in certain specific locations. 



of the outermost planet. The excess in the particle number 
with respect to /ifit(fl) is close to or above 600 particles in 
the orderly (a ~ 15.1 AU) and oligarchic (a ~ 13.06 AU) 
models, being more moderate in the other cases. We observe 
other peaks in the oligarchic growth model, at a ~ 8.35 AU 
and at a ~ 4.8 AU; in each case the excess of particles with 
respect to hfit(a) becomes more moderate as the semi-major 
axis decreases. These observations point out the existence 
of other specific locations where planets are formed pref- 
erentially in the oligarchic model, once the AMD stability 
criterion is satisfied. 

In order to show that such locations are indeed a sys- 
tematic property with full statistical meaning, in Fig. 1 10( b) 
we display the relative residuals [h(aj — ftflt(a)]/ftfit(a) of 
the corresponding fit for the four growth models. The loca- 
tion of the outermost peak in all models is at least a 65% 
off the best-fitting. Notice that the case of the oligarchic 
growth corresponds to a deviation over ~ 400% from the 
best-fitting. With regards to the secondary peaks observed 
in the oligarchic growth model, they are ~ 90% and ~ 29% 
in excess with respect to hm(a), in decreasing order of a, 
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respectively. This shows that these peaks are due to system- 
atic correlations in the data, i.e., a systematic formation of 
planets in specific locations in this oligarchic growth model; 
in other words, an enhancement of the probability to find 
planets at those specific locations. 

Figure UOf b) shows, in addition, that the peak asso- 
ciated with the outermost planet, that lies at the small- 
est semi-major axis, corresponds to the oligarchic growth 
one. This shows that migration is promoted by the com- 
bined effect of mass and velocity-selection mechanisms, as 
in the oligarchic growth model. Taking as reference the or- 
derly growth model and comparing the results for mass 
and velocity-selection mechanisms, we see that the effect is 
largely dominated by the the latter mechanism. 



4 UNIVERSALITY 

An important question in relation to the results presented 
in the previous section, is what happens if the important 
physical parameters L to t, Mdisc, and Mo have different val- 
ues. Put differently, how do the results depend upon the 
specific values of L to t, Maise an d Mo. Notice that different 
values of these physical parameters may affect the radial ex- 
tension of the disc, cf. Eq. 1101 In this section, we address 
this question and show that the results are universal. Here, 
by universality, we mean that the results of different pa- 
rameters are statistically the same under properly defined 
scaled variables. Notice that this universality may not hold 
by changing other parameters of the mod el such as e max , 
pa(a) or ft (|Hernandez-Menal I unpublished!) . This question 
is clearly of central importance for a comparison with the 
observations of exo-solar planetary systems. 

To this end, we shall define the length-scale 

-Rscalc = [£tot/Mdisc] 2 //i, (13) 

which is defined only in terms of the dynamically invari- 
ant quantities of the underlying many-body Hamiltonian 
(Eqs. [T}[2]). We emphasize here that /i = GMq includes the 
dependence on the mass of the host star Mo, and is not fixed 
to 4-7T 2 unless Mo = Mq. The dimensionless scaled variables 
are thus given by 

a = a/R Bcalc , (14) 

m = m/A/disc (15) 

According to these equations, the masses of the planets 
formed scale linearly with the mass of the disc, i.e., Mdisc 
defines simply the overall mass scale of the simulations. 
In turn, the semi-major axes scale linearly with Mo, and 
quadratically with the ratio of the total angular momentum 
and the mass of the disc. 

We compare now equivalent calculations for the same 
growth model using different physical parameters. We shall 
use the oligarchic growth model of Section 13.51 as example, 
and compare with calculations performed for the same val- 
ues of Mo and Mdisc, which therefore will not change the 
mass-scale of the formed planets, but with the total an- 
gular momentum of the system fixed to L' tot — 3.479 x 
KT 2 A/ AU 2 yr -1 . The latter value is ~ 1.57 times the value 
used in Section f3. 51 

Figure lllf a) shows the mass-semi-major axis diagram 
for this case. The figure is quite similar, in a qualitative 
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Figure 11. (a) Mass-semi-major axis diagram for the combined 
results of 4096 simulations of the oligarchic growth model with 
total angular momentum of the system set to L' tot = 3.479 X 
10 — 2 Mq AU 2 yr — 1 (Mdisc aI1( i Mq are the same as in the simu- 
lations of Section 13-56 - (b) Same as the top panel but using the 
scaled variables, Eg. 1141 and 1151 displaying the data used in (a) 
and the data used in Fig. [5Ja). 



sense, to Fig.[8ja). Quantitatively, the vertical scales of both 
figures coincide, as expected. The horizontal scale is enlarged 
by a factor close to 2.5 (actually, (Z^ ot /L to t) 2 ). In Fig jllf b) 
we display the same correlation diagram using the scaled 
variables defined in Eqs. Q3] and 1151 for the data of both 
parameter sets. The difficulty to distinguish the data from 
one case to the other clearly illustrates the meaning of uni- 
versality. We emphasize that the same scaling holds for the 
position of the peaks of the relative residuals that mark the 
systematic formation of planets at certain specific locations. 
Equivalent results can be obtained by varying the values of 
Mdisc and M . 

The universal property illustrated in Fig. [TT] is a con- 
sequence of the fact that the (N + l)-body Hamiltonian 
(Eqs. [TJQ includes only terms which build up a homoge- 
neous potential (o f degree -1), thus being in variant under ap- 
propriate scaling (lLandau fc Lifshitzlll97q '). Therefore, uni- 
versality actually holds for any dynamical variable and also 
for all the growth models discussed here. This property 
opens another perspective to statistically analyse the data 
of exo-solar planetary systems. 

Yet, the scaling property of universality may break 
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down if other parameters used in the model strongly in- 
fluence the results of the simulations. Changing the func- 
tional form of p a (a) will change the scaling properties of the 
mass in terms of the semi-major axis, which can be inferred 



from the analytical results of Laska r (2000a), cf. Eq. [9] An- 
other more subtle situation arises, for example, when a much 
larger value of L' tot is considered, keeping all other parame- 
ters fixed. In this case, the initial planetesimals disc will be 
more extended with respect to the parameter correspond- 
ing value a max , with the same functional form of the initial 
p a (a), thus making up a much fainter planetesimal disc due 
to normalization. If the density is not taken consistently, 
the value of e max used may become important in the simu- 
lations, since some collisions of particles may not take place 
for L' tot . Adjusting the number of particles at the beginning 
of the simulations yields universality again. 



5 SUMMARY AND CONCLUDING REMARKS 

In this paper, we have studied some physically-motivated 
variat ions on Laska r's simple model of accretion and evolu- 
tion (|Laskarll2000al lbl), simulating the formation and dynam- 
ics of planetary systems in the secular limit, thus excluding 
effects of mean-motion resonances. Fulfilling the AMD sta- 
bility criterion ensures that planetary collisions are no longer 
possible for the averaged system, but might be possible if 
mean-motion resonances were included. The variations we 
incorporated include important situations which are rele- 
vant to the formation of planetary systems according to the 
current underst anding (|Safronovl 1 19721 ; IPolack et al.l Il99rj; 
lArmitagdlioibl ). In particular, we implemented the model 
with a strict conservation of the total angular momentum 
of the system during the simulations. We also addressed the 
implications that arise from including in the accretion proba- 
bility distribution a dependence on the mass of the planetes- 
imals, and on the relative velocity of the colliding particles 
at the collision point. These are important factors that influ- 
ence the accretion rate processes. In our implementation of 
these physical effects we have tried to maintain the simplic- 
ity of Laskar's original model. An important limitation of 
our models, inherited from Laskar's simplified model, is the 
lack of a true dynamical evolution, which thus prevents us 
from answering questions involving or addressing the phys- 
ical time in the formation processes. 

Our statistical analysis was based on the mass-semi- 
major axis correlations of the formed systems, using the 
current values of the total angular momentum and mass of 
the planets of the Solar System as main parameters, though 
we included also remarks on the final eccentricities. For the 
comparisons, we consider as reference Laskar's model impos- 
ing the conservation of the total angular momentum of the 
system, which we referred to as orderly growth. This case 
shows a power-law behaviour of the mass of the formed plan- 
ets in terms of the semi-major axes in accordance to Laskar's 
prediction, but with a different value of the exponent. At the 
edges of the disc, our statistical results showed distinctive 
features deviating from this power-law behaviour, which are 
associated with the location of the inner- and outer-most 
planets. Introducing a mass-selection mechanisms where 
more massive bodies have an enhanced probability to ac- 
crete, thus modelling a kin of runaway growth, we found 



a comparatively smaller number of planets formed, which 
are more massive. Interestingly, we observed a systematic 
occurrence of rather small-mass remnants embedded in the 
planetary system. We also considered a velocity-selection 
mechanism by only allowing accretion when Safronov's con- 
dition is satisfied; in this case, our results showed that the 
AMD stability criterion was not fulfilled in general. We thus 
relaxed this condition, and fixed the relative orientation of 
the elliptic orbits of the colliding bodies to the value cor- 
responding to the minimum of the relative velocity. This 
criterion selects the quasi-tangential collisions, which take 
longer times of interaction, making physically plausible ac- 
cretion. In this case, we found that the average number of 
planets was further reduced, which are then more massive, 
at the end of the simulations. Furthermore, the results man- 
ifest the appearance of a gap in the mass-semi-major axis 
diagram, which defines a lower bound for the masses of the 
planets at either edge, that depends on their location. 

We also considered the combination of the mass- and 
velocity-selection mechanisms, in what we called the oli- 
garchic growth model. In this case, the number of final plan- 
ets was further reduced and their masses increased, with the 
peak reaching values over 0.6A/./. An interesting result is the 
systematic appearance of planets at specific locations, man- 
ifested as a definite excess of planets formed in this growth 
model. This property is manifested as a clustering of planets, 
i.e., an enhancement in the density profile in definite regions 
of the mass-semi-major axis diagram. We emphasize that, 
despite the stochastic nature of the model, the combined 
mass- and velocity-selection mechanisms induce strong cor- 
relations that yield this clustering. We observed from the lo- 
cation of the outermost planet, that the strongest migration 
occurs precisely in this oligarchic growth, due fundamen- 
tally to the velocity-selection mechanism, i.e., gravitational 
focusing. 

We also addressed the question of variations to the main 
physical parameters of the model. We found that the re- 
sults are statistically invariant with respect to such changes, 
whenever they are expressed in terms of appropriate scaled 
variables. The masses scale linearly with the mass of the disc, 
and the semi-major axes scale according to Eqs. [T3land[T4l 
involving the angular momentum, the mass of the disc and 
the mass of the central star. In this sense, the results are 
universal. Universality holds in particular when the same 
functional form of the initial distributions p a {a) and p e (e) 
are maintained. This property follows from the mech anical 
similarity of the Hamiltonian (jLandau fc Lifshitzlll976l ). 

In most of our simulations, we used the parameters de- 
scribing the current state of the Solar System, and an ini- 
tial homogeneous linear mass-density distribution p a (a). We 
expect that changing the initial linear mass-density distri- 
bution influences the power-law behaviour obtained, essen - 
tially as it does in the analytical results of lLaskail (|2000al ) . 
Comparing the masses of the planets formed in the simula- 
tions with those of the Solar System, in our simulations the 
inner planets are more massive, and the outer planets are 
thus not massive enough. Yet, it is interesting to remark 
that the occurrence of the specific locations (semi-major 
axes) where the planets form, obtained for the oligarchic 
growth model, are consistent with the initial conditions used 
for three (out of four) of the major planets considered in 
the simulations of the architecture of the Solar System per- 
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formed within the first Nice model bv lTsiganis et al l (|2005h : 
as mentione d, this unexpected agreement does n ot hold for 
the masses (Th ommes. Duncan fe Levisoni r2003). Consider- 
ing that our results are qualitative, such a partial quanti- 
tative agreement is encouraging to construct more realistic 
models, that fully explain the initial conditions of the Nice 
model in a statistical robust sense. 

Future work along these lines will be aimed to incor- 
porate in the model e ffects related to mean-motion reso- 
nances (lMalhotralll995l ), which are kn own to be important 
in the architecture of the So lar System (|Tsiganis et al.ll2005l ; 
iMorbidelli et al.ll2007l.l2009l 'l. Furthermore, improving on the 
assumed form of the initial linear mass-density distribution 
p a {a) may yield better comparison of the mass d istribu- 
tions. For th e concrete case of the Solar System, lLaskarl 
(1997. 2000a) suggested that a better modelling may be an 
initial linear mass-density distribution that considers sepa- 
rately the inner and outer solar systems. Other possibilities 
include the distinction of particles as gas or heavy elements, 
or the addition of an outer planetesimal disc, as assumed by 
the Nice model. 
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